Today

• Homework: post tomorrow
• More raster…
• Socioeconomic data at SEDAC
• Final Project discussion / work?

Libraries

library(knitr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(raster)
library(rasterVis)
library(scales)

Today’s question

How will future (projected) sea level rise affect Bangladesh?

1. How much area is likely to be flooded by rising sea level?
2. How many people are likely to be displaced?
3. Will sea level rise affect any major population centers?

Bangladesh

getData("ISO3")%>%
  as.data.frame%>%
  filter(NAME=="Bangladesh")

##   ISO3       NAME
## 1  BGD Bangladesh

Download Bangladesh Border

Often good idea to keep data in separate folder. You will need to edit this for your machine!

datadir="~/GoogleDrive/Work/courses/2015_UB503/SpatialR/data"

Download country border.

bgd=getData('GADM', country='BGD', level=0,path = datadir)
plot(bgd)

Topography

SRTM Elevation data with getData() as 5deg tiles.

bgdc=gCentroid(bgd)%>%coordinates()

dem1=getData("SRTM",lat=bgdc[2],lon=bgdc[1],path=datadir)
plot(dem1)
plot(bgd,add=T)

Mosaicing/Merging rasters

Download the remaining necessary tiles

dem2=getData("SRTM",lat=23.7,lon=85,path=datadir)

Use merge() to join two aligned rasters (origin, resolution, and projection). Or mosaic() combines with a function.

dem=merge(dem1,dem2)
plot(dem)
plot(bgd,add=T)

Saving/exporting rasters

Beware of massive temporary files!

inMemory(dem)

## [1] FALSE

dem@file@name

## [1] "/private/var/folders/lg/gz_jcfk5617dlpzdtg23x3rh0000gn/T/R_raster_adamw/r_tmp_2015-10-19_114157_62728_62244.grd"

file.size(sub("grd","gri",dem@file@name))*1e-6

## [1] 144

showTmpFiles()

## r_tmp_2015-10-19_105338_61539_09013 
##  r_tmp_2015-10-19_105340_61539_06688 
##  r_tmp_2015-10-19_105809_61709_12650 
##  r_tmp_2015-10-19_105934_61709_54451 
##  r_tmp_2015-10-19_105939_61709_84199 
##  r_tmp_2015-10-19_110857_61940_90130 
##  r_tmp_2015-10-19_110900_61940_66881 
##  r_tmp_2015-10-19_111013_61940_26505 
##  r_tmp_2015-10-19_111142_61940_44518 
##  r_tmp_2015-10-19_111147_61940_41991 
##  r_tmp_2015-10-19_111532_54312_02140 
##  r_tmp_2015-10-19_113231_62659_02529 
##  r_tmp_2015-10-19_113242_62673_02529 
##  r_tmp_2015-10-19_113257_62686_02529 
##  r_tmp_2015-10-19_113326_62711_02529 
##  r_tmp_2015-10-19_114157_62728_62244

rasterOptions()

## format        : raster 
## datatype      : FLT8S 
## overwrite     : FALSE 
## progress      : none 
## timer         : FALSE 
## chunksize     : 1e+07 
## maxmemory     : 1e+08 
## tmpdir        : /var/folders/lg/gz_jcfk5617dlpzdtg23x3rh0000gn/T/R_raster_adamw/ 
## tmptime       : 168 
## setfileext    : TRUE 
## tolerance     : 0.1 
## standardnames : TRUE 
## warn depracat.: TRUE 
## header        : none

Set with rasterOptions(tmpdir = "/tmp")

Saving raster to file: two options

Save while creating

dem=merge(dem1,dem2,filename=file.path(datadir,"dem.tif"))

Or after

writeRaster(dem, filename = file.path(datadir,"dem.tif"))

WriteRaster formats

rgdal does even more…

Crop to Bangladesh

dem=crop(dem,bgd,filename=file.path(datadir,"dem_bgd.tif"),overwrite=T)
plot(dem); plot(bgd,add=T)

gplot(dem,max=1e5)+geom_tile(aes(fill=value))+
  scale_fill_gradientn(
    colours=c("red","yellow","grey30","grey20","grey10"),
    trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e3)))+
  coord_equal(ylim=c(-21,25))+
  geom_path(data=fortify(bgd),
            aes(x=long,y=lat,order=order,group=group),size=.5)

## Regions defined for each Polygons

Terrain analysis (an aside)

terrain() options:

• slope
• aspect
• TPI (Topographic Position Index)
• TRI (Terrain Ruggedness Index)
• roughness
• flowdir

Use a smaller region:

reg1=crop(dem1,extent(93.8,94,21.05,21.15))
plot(reg1)

The terrain indices are according to Wilson et al. (2007), as in gdaldem.

Calculate slope

slope=terrain(reg1,opt="slope",unit="degrees")
plot(slope)

Calculate aspect

aspect=terrain(reg1,opt="aspect",unit="degrees")
plot(aspect)

TPI (Topographic Position Index)

Difference between the value of a cell and the mean value of its 8 surrounding cells.

tpi=terrain(reg1,opt="TPI")

gplot(tpi,max=1e6)+geom_tile(aes(fill=value))+
  scale_fill_gradient2(low="blue",high="red",midpoint=0)+
  coord_equal()

Negative values indicate valleys, near zero flat or mid-slope, and positive ridge and hill tops

Your Turn

• Identify all the pixels with a TPI less than -15 or greater than 15.
• Use plot() to:
  • plot elevation for this region
  • overlay the valley pixels in blue
  • overlay the ridge pixels in red

Hint: use transparent to plot a transparent pixel.

Extract peaks/ridges and valleys:

plot(reg1)
plot(tpi>15,col=c("transparent","red"),add=T,legend=F)
plot(tpi<(-12),col=c("transparent","blue"),add=T,legend=F)

TRI (Terrain Ruggedness Index)

Mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.

tri=terrain(reg1,opt="TRI")
plot(tri)

Roughness

Difference between the maximum and the minimum value of a cell and its 8 surrounding cells.

rough=terrain(reg1,opt="roughness")
plot(rough)

Hillshade (pretty…)

Compute from slope and aspect (in radians). Often used as a backdrop for another semi-transparent layer.

hs=hillShade(slope*pi/180,aspect*pi/180)

plot(hs, col=grey(0:100/100), legend=FALSE)
plot(reg1, col=terrain.colors(25, alpha=0.5), add=TRUE)

Flow Direction

Flow direction (of water), i.e. the direction of the greatest drop in elevation (or the smallest rise if all neighbors are higher).

Encoded as powers of 2 (0 to 7). The cell to the right of the focal cell ‘x’ is 1, the one below that is 2, and so on:

flowdir=terrain(reg1,opt="flowdir")

plot(flowdir)

Much more powerful hydrologic modeling in GRASS GIS